Investigation of crude oil properties impact on wettability alteration during low salinity water flooding using an improved geochemical model

In low-salinity water flooding (LSWF), modifying the injected brine composition leads to greater oil recovery from carbonate reservoirs. The processes that control improved recovery during LSWF are not totally clear, which could lead to ambiguities in finding optimum brine composition regarding wettability alteration (WA) toward water wetness. One of the methods to determine WA is bound product sum (BPS) calculation using geochemical tools. In the case of wettability improvement, the BPS value of a crude oil-brine-rock (COBR) system should be at its minimum value. In this study, an improved geochemical model is developed, which includes the effects of oil composition (i.e., acid number, base number, and weight percent of nonhydrocarbon components) and physical properties of oil (i.e., density, viscosity, and solution gas-oil ratio) on COBR interactions. The proposed method generates BPS as a function of temperature, pressure, oil and brine composition, and pH for carbonate rocks. The model applicability was validated using several experimental data sets available in the literature. The results of the improved BPS model were in line with the results of contact angle and zeta potential measurements as the major indices of rock wettability. BPS calculations using the available geochemical tools sometimes failed to predict the correct WA trend since they overlooked the impact of oil properties on COBR interactions. The model predictability was also compared with the results of an available geochemical tool, PHREEQC, and the results demonstrate just how important the effect of oil properties and composition inclusion on wettability determination is. The improved BPS approach could be successfully utilized as an optimization tool to optimize the water composition during LSWF for a given COBR system.

Phase equilibrium. The equilibrium between the aqueous phase and oil phase is evaluated by equating the fugacity of each component in the two phases. The fugacity of species in the liquid phase is calculated using the Peng-Robinson equation of state assuming an initial composition 30 . The fugacity coefficient yields the fugacity using oil composition and pressure, The fugacity of the species in the aqueous phase is calculated using Henry's law 31 : www.nature.com/scientificreports/ where H i is Henry's law constant for component "i" at pressure P and temperature T. In this equation, H s i is Henry's law constant at water saturation pressure P s H2O , is the gas constant, and υ i . is the partial molar volume of component i in the aqueous phase at T 31 .
Thus, the objective function is written as follows: where i refers to common species that can transfer between phases (in this model, N 2 and CO 2 are considered soluble components in the oil and water phases).
Aqueous species reactions. The equilibrium state in the aqueous phase can be calculated by equating the equilibrium constant of the reaction and the activity product of the species: where α is the number of aqueous reactions, Q is the activity product of the reaction (aA + bB = cC + dD), in which A and B are reactants, C and D are products, a, b, c, and d are stoichiometric coefficients in the reaction, and ac is the activity coefficient of the aqueous species in the reaction:

Surface complexation reactions
Rock/brine surface Oil/brine surface www.nature.com/scientificreports/ In this model, the activity could alternatively be calculated using a modified Debye-Huckel (also termed the Bdot model, up to an ionic strength of 0.7 mol/kg) or Pitzer model 32 . The equilibrium constant calculation is based on the enthalpy difference, volume difference of reaction, and analytical model which is used in PHREEQC software 28 . Mineral reactions. The reaction rate of minerals is calculated through the following equation 33,34 : where m indicates mineral reactions, Q is the activity product of thmineral equation (the activity of minerals is assumed to be unity) and A β is the specific surface area of the rock. The equilibrium constant of mineral reactions (K eq,m ) is calculated using the analytical expression used in PHREEQC and using the PHREEQC database 28 .

Surface complexation reactions.
In the surface complexation model, the calcite surface consists of chemically active surface sites that react with ions of the aqueous phase(electrolyte). The strength of these interactions is described by the equilibrium constants. The equilibrium concentration of surface complexes shows the net surface charge and the surface potential 35 . In the rock/brine surface, two main reaction sites are considered > Ca and > CO 3 , where the > sign indicates complexes at the rock surface. Additionally, in the oil/brine contact, two main sites are considered -NH and -COO, where the-sign indicates complexes at the oil surface. The fraction of each surface complexation at each surface site should be calculated using the equality of the equilibrium constant and activity product of surface reactions by considering electrostatic interactions on surfaces 36 . For instance, in the case of the rock/brine surface complexation reaction of R-16 in Table 1, the equilibrium constant is calculated as follows: where g is the mole fraction of surface complexation species, e.g., for species "i" at the rock or oil surface 36 : where N i is the number of moles of complexation species formed at the rock or oil surface. ρ is the bulk rock density or oil density, which is calculated at the desired pressure and temperature employing empirical correlation using the properties of stock tank oil, and s d is the surface site density of rock or oil (i.e., number of surface sites per unit area of the rock surface), which has been estimated to be between 2 and 8 sites/nm 2 for carbonate rocks 37 . The site density of oil depends on several parameters, such as acid and base number, and A β is the specific surface area of oil or rock. The surface area of oil is proposed to be 0.1 m 2 /g 27 or 0.11 m 2 /g 19 . The cations of high charge density (i.e., smaller and highly charged) are more hydrated and reluctant to bind to the rock surface and this has a significant role in wettability alteration 38 . The surface area of oil/brine depends on some parameters, such as the acid number/base number ratio and brine composition, because in carbonate rocks, the reactivity of carboxylate groups (-COO) is more important than that of -NH groups since N does not form a stable complex with aqueous species. The residuals of the surface complexation reactions are as follows 39 : The electrostatic interaction term (χ) in the equilibrium constant of surface complexation reactions is: where F is the Faraday constant, R is the universal gas constant, is the net charge over the surface complexation reaction and ψ 0 is the surface potential. The surface and zeta potentials at a distance from the surface in the diffuse layer are related by the Gouy-Chapman theory 40 .
The magnitude of the Debye length (a point where ions in the bulk electrolyte solution do not feel electric attraction, the thickness of this diffuse layer is called the Debye length) depends solely on the properties of the solution, not on any properties of the surface, sucas its charge or potential 41 .
It is noticeable that the molality of species is based on the weight of water in the aqueous section of the thin layer between oil and rock. Hence, the portion of water that is involved in the diffuse double layer must be excluded in equilibrium calculations. By defining cation exchange capacity 36 , the mole fraction equation is obtained as follows: For rock/brine surface species, cation exchange capacity, related to the mole fraction of surface species by the above equation, is considered for surface charge calculation as 42 : where Sw is the saturation of water at the given grid block. The enrichment factor is useful for modeling the relative enrichme or depletion of equally charged species in the electrostatic layer on a charged surface, which is related to enhanced complexation in a low dielectric permittivity medium 43 .
The acid number and base number strongly affect the oil/brine surface compositions 44 . The base number to acid number ratio (BNANR) is introduced as an enrichment factor in the model. The base number represents the -NH group and the acid number represents the -COOH group. Therefore, the presence of -NH to -COOH groups at the oil/brine surface is related to BNANR.
For the rock/brine surface species, the phase proportion is usually considered to be 0.1 for simulation applications 32 . This means that 0.1 of the moles of minerals participates in the formation of surface complexations.
The equation set includes the mentioned chemical and thermodynamic equilibria, volume constraint equation (i.e., the total volume of the phases should be equal to the pore volume of the rock), and mass balance equation of the species (i.e., the mass of any species in the pore volume should be equal to the rate of change of its mass). More details on the mass balance and the volume constraint equations are available in the work of Nghiem et al. 45 . The root of the set of nonlinear algebraic equations representing the whole geochemical system at equilibrium should be solved. The system of equations is solved using a trust-region-reflective algorithm and gives the equilibrium concentrations of species. A summary of reactive transport code functions is illustrated in Fig. 1. The application of the functions is presented as follows: Main is for the first guess calculation, advancing time steps, and plotting the results, t0Calc is for first guess calculation using initial values in Initialguess0 function, Initialguess2 is for the calculation of iterative time step values, Activity calculates activity using Bdot or Pitzer model, ODEN calculates oil density using oil API and dissolved gas-oil ratio, Viscosity calculates the viscosity of brine using salinity (or TDS) and oil at the desired pressure and temperature using physical properties of oil (i.e., oil density, dissolved gas-oil ratio, and bubble point pressure), PhaseEql is a function of phase equilibrium calculations using Peng-Robinson and Henrys' constant calculations, corey returns the relative permeability of oil and water and capillary pressure using Brooks-Corey's equilibrium model, calc_logk calculates the equilibrium constants of reactions through the analytical model of PHREEQC, ddl calculates the portion of water that exists in the diffuse double layer and must be excluded in the concentration calculations based on the mass of water.
The wettability of the COBR system is determined using a surface complexation model that is considered an improved BPS approach. The flowchart of the developed model for wettability determination using BPS calculation is presented in Fig. 2. In the following section, the BPS calculation of a given COBR system is discussed.

Surface complexation model
Surface complexation modeling presumes an electrical double layer at each interface and the existence of charged surface species whose concentrations depend upon the chemical makeup of water, oil, and mineral surfaces 33,46 .

Oil-brine surface charge
The diffuse layer model is a simple model for the electric double layer and envisions a surface layer with a pHdependent surface charge and a Gouy-Chapman diffuse layer of oppositely charged ions; see, e.g., Dzombak and Morel 47 , with the sorbed ions residing on the surface layer.
Measured zeta potentials of the oil-brine interface could be represented through the sum of carboxylic acids and nitrogen bases 26 . In the oil/brine surface reactions, -COOH and -NH represent dangling carboxylate and nitrogen base groups present at the oil-water interface. At pH > 4.9, any rise in the ionic strength will increase the number of deprotonated carboxylate groups and interfacial negative charges, while at pH < 4.9, any rise www.nature.com/scientificreports/ in the ionic strength will increase the number of protonated nitrogen bases and positive surface charges. The absolute value of oil-brine interface zeta potential is decreased by increasing the ionic strength of the solution 27 . Relative numbers of acid and base groups in the model may be changed to consider oils with different acid and base numbers. These parameters, along with many others, will affect surface complex concentrations involving reactions of a COBR system that controls the wettability of rock. LSWF is successful only if the zeta potential changes at the mineral-brine interface become more negative (i.e. decrease ionic strength) when the oil-brine interface has a negative zeta potential in the formation brine 48 . A quantitative measure of electrostatic attraction is termed BPS 27,32 . This electrostatic attraction causes oppositely charged species to attract each other, leading to an oil-wet surface. BPS is a reliable index that indicates the measure of the electrostatic bond between oppositely charged surface complexes, which could also be an indicator of wettability or tendency of oil to stick on the rock surface; therefore, BPS contains electrostatic charge and surface potential and interaction of different species together.
The lower limit of BPS is "zero", meaning there is no electrostatic bound between surface species. Therefore, low values of BPS indicate a water-wet system of COBR. In this model, BPS is calculated as follows (absolute unit is mol 2 and density unit is (μmol/m 2 ) 2 ): BPS values should decrease for wettability improvement (i.e., more water-wet state), based on Brady et al. 26 . A schematic of the electrostatic interaction between ions in a COBR system is shown in Fig. 3.
The system of nonlinear equations is solved using a trust-region-reflective algorithm, so the equilibrium state of reactions is determined. First, it should be clarified that the wettability alteration trend of a given COBR system is predicted correctly by employing the developed model, so the model validation is performed using several experimental data sets in the literature that are discussed in the following section. Next, method validation is Calculate bond product sum of COBR system( BPS)

Calculate equilibrium concentraion of species( includes aqueous, mineral and surface species)
Solve system of equations

Model validation
Two experimental studies were used to validate the results from the model.

Effect of injected brine composition on calculated BPS. The contact angle is a measure of wettabil-
ity, which is the angle of the fluid-fluid interface. If the contact angle of one fluid is smaller than 90°, the solid is said to be preferentially wet by this fluid 49 . Alshakhs and Kovscek 18 conducted zeta potential measurements on calcite rock (different pH values of injected brine) to calculate contact angle employing DLVO theory using disjoining pressure calculations. The petrophysical properties of the core samples used for core flooding are shown in Table 2, and the properties of the crude oil are shown in Table 3. The ionic composition of different injected brines is shown in Table 4. The temperature of the experiments was 60 °C.
At pH = 8.5, the measured and calculated contact angle is shown in Table 5. Sequence of contact angle values of injection scenarios at this pH is: The estimated BPS[mol 2 ] values using the model are shown in Table 6. The sequence of BPS values of injection scenarios at this pH is BPS 1 > BPS 4 > BPS 3 > BPS 2 , which implies that BPS correctly models the wettability alteration during injection of different composition brines and indicates that the presence of MgSO 4 2− in the composition of injected brine has a positive effect in altering wetness toward water wetting, while SO 4 2− has very little effect. In the case of injection scenario 4 relative to 1, it is clear that    www.nature.com/scientificreports/ negatively charged SO 4 2− attracts the positively charged calcite rock surface; therefore, there are more water-wet sites for the negatively charged oil surface. Injection scenario 3 causes calcite dissolution due to the higher activity of Mg 2+ relative to rock Ca 2+ ions, and consequently, injection scenario 2 results in more wettability improvement due to the simultaneous effect of SO 4 2− attraction and dissolution of calcite. It seems that the improved model captures the mentioned effects concurrently since there is a satisfactory agreement between the model results and the experimental results. The calculated BPS values using the model are presented in Fig. 4. By decreasing PH of the solution, concentration of H + ions increases, causing the reaction between divalent ions of formation brine (Ca2 + , Mg2 +) with the solution to decrease 50 . This results in concentration decreasing of carboxylate groups (-COO − ) on oil surface which leads to water-wet state. Consequently, by consideration of positively charged calcite surface, the sum of oppositely charged complexes at the oil-brine-rock interfaces decreases that means BPS reduction.
Effect of temperature on calculated BPS. The developed model was also validated using the results of Mahani et al. 19 experiments. They conducted zeta potential and contact angle measurements on several carbonate rock samples at different temperatures. The injected brine was 25 times diluted seawater, and its composition is shown in Table 7.
The calcite site density and specific surface area are considered 2 sites/nm 2 and 0.11 m 2 /g, respectively. The physical and chemical properties of crude oil A are shown in Table 8.
It was reported that in limestone A (only calcite rock), increasing the temperature did not show any wettability improvement, whereas dolomite (99% dolomite + 1% quartz) showed wettability improvement by injecting 25 times diluted seawater (25dSW).
The calculated BPS values using the developed code are presented in Table 9. www.nature.com/scientificreports/ As shown in Table 9, in the case of dolomite rock type, temperature rise causes improvement in wettability through water wetting, whereas in the case of calcite rock, increasing the temperature has a negligible effect on wettability. The calculated BPS values of Mahani et al. 's work 19 and the average contact angle changes of the two rock samples (limestone A and dolomite), reported in his work, are illustrated in Fig. 5. The contact angle change is the difference between the contact angle of an oil droplet at the time of equilibration with formation water and the end of exposure to 25dSW.
Different mechanisms for wettability alteration have been proposed in the literature 4,36,51 that include altering the surface charge through adsorption of Ca 2+ and SO 4 2− on the chalk surface and Ca 2+ substitution by Mg 2+ due to increasing ion reactivity at higher temperatures. By injecting brines containing SO 4 2− , it is assumed that this ion will be adsorbed on a positively charged chalk surface. Thus, the bond between the negatively charged oil surface and the positively charged surface will break down. This effect represents a correlation between oil recovery and temperature in spontaneous imbibition tests, through intensifying this effect by increasing temperature 52 . It is proposed that at high temperatures, the reactivity of these ions increases, and Ca 2+ on the rock surface will be substituted by Mg 2+ .
These observations are thoroughly included in modeled reactions through activity coefficient calculations, equilibrium constant databases, and molar volume changes in the reactions, which are functions of pressure and temperature. However, to the best of the authors' knowledge, the effect of increasing pressure on wettability alteration has not been presented in the literature before use in model validation. The effect of pressure on the reaction is incorporated through the molar volume change of the reaction that affects the equilibrium constant calculation, although the fugacity of the phases is affected by pressure. Therefore, the pressure effect is also incorporated into the phase equilibrium calculations. This finding confirms the usefulness of the new model as an optimization tool that optimizes injected brine composition by considering the effect of various parameters, such as pressure, temperature, and oil properties. On the equilibrium concentration of all species.
Measured zeta potentials, due to the evaporation of brine at high temperatures, are much more stable at 25 °C. Hence, zeta potentials at 25 °C are chosen to examine the surface potential dependencies of oil and rock surface reactivity at higher temperatures, assuming that zeta potentials tend to decrease uniformly with an increase in temperature 33 . There are a few measured zeta potential data in the range of oil reservoir temperatures (higher than 100 °C), and the existence of a model could be of assistance to optimize brine compositions and be employed to inject low salinity water flooding under reservoir conditions.   51 crude oil is 1 mg KOH/g of oil, which was not reported by the authors. The temperature of the experiments was 100 °C. All other properties regarding core, oil, and brine were the same as those reported in their work. For the examined dilutions of seawater presented in Yousef et al. 51 , BPS was calculated in density units, and the results are shown in Fig. 6. The reported contact angle by Yousef et al. 51 is shown in Table 10.   www.nature.com/scientificreports/ It is clear that since BPS has a little positive response (i.e., wettability improvement or BPS reduction) compared to formation brine, injecting seawater has little effect on wettability. However, two-and tenfold dilutions of seawater (2dsw and 10dsw, respectively) have a significant effect on the BPS trend since at a constant pH, seawater causes smaller BPS values than formation brine. It is also reported that 20 times and 100 times dilutions of seawater (20dsw and 100dsw, respectively) have a lesser effect on the BPS trend in comparison to the previous dilutions (10dSW and 20dSW, respectively). These observations correspond to the experimental work by Yousef et al. 51 .
Dilution of seawater leads to ionic strength reduction, which is directly related to the ionic activity coefficient. Thus, the interaction between ions decreases, which causes interfacial tension reduction and subsequent wettability improvement. As already mentioned, in the existing geochemical tools, such as PHREEQC 28 , the effect of oil phase properties and compositions is excluded, whereas, in the developed model, it is included in the equilibrium calculations. This comparison emphasizes the method validation and indicates that the presented model could effectively determine wettability alteration relative to the existing geochemical tools.

Conclusions
This study presents an improved BPS approach for determining wettability alterations of carbonate rocks by calculating the BPS values of a given COBR system. The effect of oil properties (i.e., density, viscosity, and solution gas-oil ratio) and composition (i.e., acid number, base number, and weight percent of nonhydrocarbon components) is included in COBR interactions. Reactions between phases and all aqueous, mineral, oil, and surface species are coupled in the developed geochemical model, and the system of equations is solved to give equilibrium concentrations of species. BPS values are calculated using the computed equilibrium concentration of species for a given COBR system. The results of the developed model are validated using several published experimental works. The following conclusions were drawn from the study: 1. The developed program considers the effect of oil composition and its physical properties on COBR interactions and wettability. 2. The proposed method conducted phase equilibrium calculations, while it considered the effect of soluble components of each phase on the composition and thickness of the diffused double layer (since the portion of water that is involved in the diffuse double layer must be excluded in equilibrium calculations). Subsequently, it generates BPS as a function of temperature, pressure, oil and brine composition, and pH for carbonate rocks. 3. The developed model is validated using experimental data in the literature; contact angle alteration in the two-phase oil-brine system of different brine compositions showed the same trend as experimental results using zeta potential measurements. 4. Contact angle alteration in the two-phase oil-brine system at different temperatures showed the same trend as the experimental results using zeta potential and contact angle measurements. Therefore, the developed model could successfully determine wettability alteration by considering the temperature effect on the equilibration of phases and species. 5. The results of BPS calculations using PHREEQC are compared to the results of the developed program for the oil-brine-rock system of experimental work by Yousef et al. 24 . This result supports the idea that oil phase properties and composition inclusion are of great importance in phase and species equilibrium calculations. 6. The results suggest that the improved BPS approach could be successfully utilized as an optimization tool to optimize the water composition during LSWF since BPS (which indicates the wettability) is calculated as a function of pressure, temperature, oil and brine compositions, and pH for a given COBR system.

Data availability
The data will be available upon request.